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MOLECULAR DYNAMICS OF POLAR GAS ROTATIONAL RELAXATION 
by Frank! Zeleznik, John V. Dugan, Jr., and Raymond W. Palmer 

Lewis Research Center 

SUMMARY 

A 50 -molecule sample in contact with a heat bath of similar molecules was used as a 
model to calculate the rotational collision numbers for HCl and DCl at 300 and 

500 K, by the method of molecular dynamics. The mean rotational energy of the sample 
exhibits an exponential approach to a steady -state value; however, there was considerable 
scatter about the exponential curve. The results for generally agree with previous 

theoretical calculations and these results are intermediate to both experimental acoustic 
data and the earlier perturbation results. 


INTRODUCTION 

Accxirate experimental studies of relaxation phenomena on the macroscopic level are 
difficult at best. This difficulty is compounded by the existence of two different theoret- 
ical definitions of a relaxation time and, furthermore, these definitions are generally 
used interchangeably. Previous papers (refs. 1 and 2) applied these two definitions to 
the special case of rotational relaxation in polar gases. Both of these papers were essen- 
tially analytical calculations and were motivated by the desire to obtain a closed-form ex- 
pression for the rotational collision number Z^.^^ arising from each definition. In both 
cases a classical perturbation calculation was carried out through third order. Refer- 
ence 1 used the energy relaxation definition of a relaxation time, while reference 2 em- 
ployed the definition based on volvune viscosity. 

The results of references 1 and 2 can be succinctly stated in the following manner: 

(1) The two definitions agree exactly through third order in perturbation (apart from 
a multiplicative adjustable parameter). 

(2) Theoretically calculated values of increase for either a decrease in the 

dipole moment or a decrease in the moment of inertia; generally speaking, Z^^^ in- 
creases with temperature, but in some cases it may exhibit a shallow minimum. 
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(3) The extraction of experimental rotational relaxation times from polar gas ther- 
mal conductivity measurements is highly suspect. 

These conclusions can only be given tentative acceptance since they are based on a cal- 
cxilation which was not only a perturbation calculation but also was restricted to a two- 
dimensional model. 

The present report examines the effect of both the perturbation calculation and the 
two-dimensional model vis-a-vis the energy relaxation definition of a relaxation time. 

All the results are for a three-dimensional model of a polar gas whose rigid -rod mol- 
ecules interact by means of spherical hard cores with embedded point dipoles aligned 
along the symmetry axis. Further, all calculations are carried out numerically to avoid 
the perturbation assumption. The molecular dynamics approach to studying rotational 
relaxation in polar gases is used. 

Generally, the method of molecular d 3 mamics is applied to dense fluids or solids to 
study both their time -dependent and time -independent properties. For example, Rahman 
(refs. 3 and 4) studied liquid argon by this technique and calculated the pair -correlation 
function, the self -diffusion coefficient, the velocity auto -correlation function, and the 
time -dependent pair -correlation function. Similarly, Verlet (refs. 5 and 6) calculated 
the temperature, pressure, internal energy, and the pair -correlation function and its 
Fourier transform for liquid argon. In these cases the system "thermalized" qxiite rap- 
idly and required on the order of 2x10 seconds to reach equilibrium (refs. 5 and 7). 
Kohler, Bellemans, and Gancberg (refs. 8 and 9) studied a system of weakly interacting 
electric dipoles confined to a rigid, square lattice. They examined the approach to equi- 
librium for such a system in terms of the first, second, and fovirth moments of angular 
momentum and Boltzmann's H -function - the second moment being essentially the kinetic 
energy. They experienced difficulties in extracting relaxation times because of large 
fluctuations. An interesting aspect of this work is the demonstration of the existence of 
a set of initial conditions for which the system does not exhibit a definite temporal evolu- 
tion. This is similar to the situations observed by Hirooka and Saito (ref. 10) in their 
study of a two-dimensional square lattice of coupled anharmonic oscillators. They dem- 
onstrated the existence of an induction period that preceded the evolution toward equilib- 
rium. For a sufficiently weak coupling of the oscillators they could not detect any appre- 
ciable approach to equilibrium. 

All the examples of the molecular dynamics approach discussed to this point were 
similar in the sense that a macroscopic system was approximated by a finite system. 

The classical equations of motion for this finite system were then solved numerically. 
Such a procedure is necessary for liquids and solids because of the high densities in- 
volved. For low -density gases, on the other hand, the binary collision approximation 
can be invoked. This permits the decoupling of the large set of equations describing the 
evolution of the model into smaller subsets, with each subset now describing the interac- 
tion of a pair of molecules. It is possible to simplify the procedure even further by 
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introducing the idea of a heat bath, an infinitely large system that always remains in 
equilibrium. If a small sample system is immersed in the heat bath, all interactions 
among members of the small system may be ignored and only interactions between the 
heat bath and the sample system need be considered. A model of this kind has been used 
to study the relaxation of diatomic molecules in a heat bath by Alterman and Wilson (ref. 
11), Raff (ref. 12), and Berend and Benson (ref. 13). 

Alterman and Wilson studied the vibrational relaxation of 50 bromine molecules in a 
heat bath of xenon atoms. They terminated their calculations after fewer than 250 co- 
linear collisions between the molecules and atoms. A vibrational relaxation time was 
obtained from the initial slope of the mean vibrational energy of the 50 molecules as a 
function of the number of collisions. It was not possible for them to determine relaxation 
times from all the calculations since in some cases there was no obvious choice of an 
initial slope. 

The rotational relaxation calculations of Raff (ref. 12) and Berend and Benson (ref. 

13) are similar in gross features; however, they differ considerably in their fine struc- 
ture. Both calculated a relaxation time in terms of a mean rate of energy transfer which 
effectively plays the same role as the initial slope in the Alterman and Wilson (ref. 11) 
calculations. However, Raff used a three-dimensional average calculated over a three- 
dimensional flux distribution, while Berend and Benson performed a two-dimensional cal- 
culation and averaged over a two-dimensional Boltzmann distribution. In view of these 
differences it is not surprising to find that their results do not agree when both calculate 
a common system. Thus, Raff calculations are lower by a factor of 5, while Berend 
and Benson fe calculations are higher by about the same amoimt when compared to the 
acoustic data of Jonkman and Ertas (ref. 14) for the rotational relaxation of P-H 2 by He. 

The calculations described in this report are more akin to those of Alterman and 
Wilson (ref. 11) than they are to those of Raff (ref. 12) and Berend and Benson (ref. 13), 
in the sense that we follow the relaxation of a polar gas coupled to a heat bath. However, 
in contrast to the one -dimensional treatment of Alterman and Wilson, we shall carry out 
a three-dimensional calculation for a sufficient number of collisions to assess whether or 
not there is an exponential approach to eqmlibrium. Also we shall try to extract rota- 
tional collision numbers as well as their dependence on the temperature, dipole moment, 
and moment of inertia. Specifically, a relatively small sample containing N molecules 
is placed into a heat bath of the same kind of molecules. An observation of the rotational 
energy of the sample as a function of the number of collisions n between the sample and 
heat bath molecules should enable us to determine the rate at which the sample ap- 
proaches an equilibrium with the heat bath. If the sample approaches equilibrium expo- 
nentially, we should be able to extract a relaxation time from the data, provided that the 
sample and the heat bath are always in translational equilibrium but are not initially in 
rotational equilibrium. The rotational relaxation time for such a situation is defined 
by 
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( 1 ) 


dER_ 

dt Ntjj 

where Ej^ is the rotational energy per particle of the sample and Ej^(t = oo) = E^. If 
ER(t = 0) = Eq , integration of equation (1) gives 

/ -t/TpN\ 

Er - = <^00 - Eq) - e ^ J (2) 

If the mean time between collisions is denoted by t. , then the number of collisions n 

c 

and the rotational collision number are defined by 


n = 


( 3 ) 


7 

“^rot — 
^c 


Combining equations (2) and (3) produces 

Er - Eq = (Eoo - Eo) 


1 - exp ( - 


n 


NZ 


rot/J 


( 4 ) 


Expression (4) can be used to extract rotational collision numbers from the variations of 
Er with n. 


ANALYSIS OF THE CALCULATION 

The numerical experiment which we shall carry out is conceptually quite simple, but 
its interpretation is far from simple. To see this, consider an ensemble of such exper- 
iments; that is, we construct a very large number of identical replicas of the N- 
molecule sample system. Each of these replicas is observed as its rotational energy 
evolves when it is placed into contact with a heat bath for a fixed time period (or a fixed 
number of collisions). Because of the random natvire of the interaction between the sam- 
ple and the heat bath and because of the finite sample size, it is unrealistic to expect 
either that the evolution of Er will be a smooth exponential or that all members of the 
ensemble will evolve along a common curve. Thus our ensemble of experiments gener- 
ates an ensemble of curves and each cimve exhibits many fluctuations or changes in slope. 
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No one member of the ensemble of numerically generated ’’experimental*’ curves 
adequately describes the real physical situation; however, it is reasonable to expect that 
an ensemble average of the data would be physically relevant and could be used to extract 
a (see top half of fig. 1). Further, the fluctuations in the data for one member of 

the ensemble about the mean of the ensemble data are minimized by a sufficiently large 
sample size N. Hence, the ideal, but economically unfeasible, numerical experiment 
would be to observe an ensemble of samples, each sample consisting of very many mol- 
ecules, for a very long period of time. By contrast, the practical calculation would be 
to observe one member of an ensemble of samples, each sample consisting of a few 
molecules , for a relatively short period of time. Considerations of computation time 
alone require that N be about 50, certainly not much larger. The length of time that a 
sample must be observed in order to adequately characterize its temporal evolution is 
certainly open to question. Obviously, the observational period must be considerably 
longer than the periods of the fluctuations in the curve. Intuitively, it seems reasonable 
to reqxiire that the observational period at least satisfy the inequality n > 

Let us now assume that we have observed, for a relatively short period, the evolu- 
tion of Ej^ for one sample of size N. How is this data to be analyzed, in terms of equa- 
tion (4), in order to obtain a physically relevant Z^^^? We can either regard equation (4) 
as containing two fitting parameters Z^^^ and E^, or else we can fix E^ at the value 
corresponding to the equilibrivim rotational energy E^^ and use only as a fitting 

parameter. Presiunably, the two methods would give the same result for Z^^^ when ap- 
plied to the ensemble data of the ideal experiment discussed previously. It is, however, 
unlikely that they will give the same results for Z^^^ when applied to the practical cal- 
culation. It is not possible to say which method will give the more physically realistic 
result, although the two -parameter method will certainly give the best possible represen- 
tation of the temporal evolution of Ej^. For this reason both the two -parameter and one- 
parameter methods are used to analyze the data. Within either the two-parameter or 
one-parameter framework we must still attempt to obtain a Z^^^ from the data of one 
member of the ensemble which might be comparable to that determined by the data of the 
entire ensemble. For this reason a heuristic alternative to the analysis of the ensemble 
data is described. The description is given in terms of the two -parameter treatment of 
the data, but it is equally applicable to the one-parameter treatment. 

The random nature of the interaction between the heat bath and the sample has been 
previously pointed out. Thus, if equation (4) were fitted to the data of each ensemble 
member separately, a somewhat different pair of ”best” parameters would be obtained, 
reflecting the different history of each replica (see bottom half of fig. 1). In effect, 
the parameters Z^^^ and E^^, which characterize the evolution of the replica during 
the short observational period, are random variables. In principle, it should be possi- 
ble to construct new random variables from Z^.^^ and E^ which will serve as unbiased 
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estimators of the physically relevant and which would be obtained from an 

analysis of the data of the entire ensemble. 

Now observe that the duration of observation simply cannot remove the randomness 
of the heat bath - sample interaction. Suppose then that a given member of the ensemble 
has been observed for a sufficiently large number of collisions n^^ > and the data 

analyzed to extract the two parameters and E^(l). This sample is allowed to 

continue its evolution until there have been an additional An collisions and then the 
parameters Zj,^^(2) and E^(2) are obtained. Proceeding in this fashion, Z^^^(k) and 
E^(k) are obtained at n^^ = + An collisions. Each fixed increment An is a random 

increment in the random history of the sample and, therefore, it appears plausible to ap- 
proximate the ensemble average Zrot by the mean of ^rotW 

/ \ 1 i”? 

r>-°Vphy3ical “ ^ ^ 

mdiA ^ 

The value of k is determined largely by what is considered to be an economically 
acceptable amount of computation. 

The mean of the values of E^(k) can be used in an attempt to justify the termination 
of the calculation, a posteriori, in the following manner: It is well known that the mean 
of a random sample of size N drawn from some distribution is not, in general, equal to 
the first moment of the distribution. Instead, a set of such means is distributed about 
the first moment. It might be expected that the mean values of E^ from calculations 
done at different conditions would exhibit variations about E^^ comparable to those ob- 
tained in drawing a random sample of size N from an equilibrium distribution. On this 
basis, then, if the distribution of the mean values of E^ is radically different from that 
expected of random samples of size N, we would suspect that an insufficient number of 
collisions had been studied. On the other hand, we could reasonably hope that if the dis- 
tribution of mean values of E^ is consistent with that expected of a sample of size N, 
then enough collisions had been studied. For this purpose we must know the probability 
that a sample of size N, when drawn from an equilibrium distribution, will have a rota- 
tional energy within a given interval about the population mean. 

Consider the equilibrium distribution of molecvilar rotational energy e. If the mol- 
ecule possesses two degrees of freedom, the distribution function at a temperature T is 
given by 


f(€)d€ = iSe"^^ de 0 ^ e (5) 

-1 9 

where ^ = kT and k is Boltzmann constant. This is an example of a x - 

distribution with v degrees of freedom 
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f(X^)d(x^) = -i- ^ «-xV2h^^2 


2^^ 


1 


d(x") 


( 6 ) 


where v has been set equal to 2. This can be easily seen by making a change of varia- 
bles = 2^e in equation (5). The moment -generating fvinction for equation (6) is given 
by 


M 2(e) 


•'f\ 


f(x^)d(x^) = (1 - 2e)“*^/2 


(7) 


Instead of considering the distribution of 


N 


E 


R 






( 8 ) 


it is convenient to consider the distribution of the dimensionless quantity /3Ej^. 

N 


/3Ej,= 


- E ^^i 


(9) 


From the properties of the moment -generating function (ref. 15, p. 135), it follows that 




/■ 

*/0 


/3Ejje 

e ^ F(/3Ej^)d(/3Ej^) = 


M 


/8e| 


,N 


( 10 ) 


where F is the distribution fvmction for /3Ej^. But using the transformation x^ = 2/3e, 
as before, M^^(0) can be expressed in terms of M 2(e). 




( 11 ) 


Hence, combining equations (7) and (11) with equation (10), 
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-Nv/2 


where v should actually be set equal to 2 Thus, /3Ej^ is also distributed as a x^- 
distribution but with Nv/2 degrees of freedom In fact, it can be verified that /3Ej 
has the distribution 

F0Ejj)dCi8Ejj) = exp (- N^Ej^)d(i3Ej^) 


The moments of this distribution are 


-n /NiA/Ni' 




2 A 2 


+ n - 1 


while its cumulants are 


In (e) = nI-"" (n - 1): (15) 

O"" " e=o " 

Note that, for large N, /3Eg is a normal variable, with = ■^lJv/2. 

Now the probability P that a random' sample of size N will have a /3Eg that lies 
in the interval (1 ± p)/^^ about the mean can easily be calculated. Here p satisfies 
the condition 0 < p < 1. 


»(l+p)|Ui 




F(^Eg)d(^Eg) = 




^Ni<l+p)/2 

I 

Ml 

\ 2 /•'NKl-p)/2 


,(N^/2)-le-^ dz 


The probability P can be readily evaluated in terms of the incomplete F-function ratio 


I(u,q) H 


/.uyq+1 

my / 

•'0 


v^e"'^ dv 
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which has been tabulated by Harter (ref 16). 



Figure 2 is a plot of P as a function of Nt/2, with p taken as a parameter. It is ap- 
parent from this figure that the one-parameter family of curves grows quite slowly be- 
yond Ny/2 = 50. Since v equals 2 for the rotational energy distribution, the sample 
size N would have to be quite large for the sample means to be sharply clustered about 
the first moment. Since N has been chosen equal to 50, there should be substantial 
fluctuations in the mean values of E^. 

CALCULATION METHOD 

The calcTilational method has been described in general terms. Some of the specif- 
ics are presented in this section. The physical system is a gas composed of polar mol- 
ecules which are represented as rigid rods surrounded by a spherical hard cores with 
embedded point dipoles. The rotational energies of a sample of N = 50 such molecules 

O 

is distributed randomly over a x -distribution with v = 2 degrees of freedom and a tem- 
perature of T = Tq. The translational velocities of these molecules are assumed to be 
in equilibrium with a heat bath of the same molecules and are characterized by a 
Boltzmann distribution at a temperature T > Tq. Thus, it is not necessary to keep 
track of the translational energies of the sample molecules diiring the computation. The 
calculation is carried out by performing the following steps repetitively: 

(1) One ’’sample” molecule and one ’’heat bath” molecule are randomly selected as 
collision partners. 

(2) The two molecules scatter and, in so doing, alter the rotational energy of the 
’’sample molecule. ” 

(3) After each scattering, a new mean rotational energy of the sample molecules is 
computed. 

The details of random selection are given in the following section; the balance of this sec- 
tion is devoted to a discussion of the scattering equations of motion and their numerical 
solution. 

The relative motion of two molecules is described by the Lagrangian equations ob- 
tained from the Lagrangian 

L = Ttrans ^ ‘ V - 
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where is the kinetic energy of relative motion, Tj.Q|.(i) (where i = 1,2) is the 

rotational kinetic energy of the i^^ molecule, V is the dipole-dipole interaction, and 
is a spherical hard-core potential of diameter <r = (Uj + o^/2 (where o ^ and Og 
are the molecular diameters of the colliding molecules). The rotational kinetic energy 
is conventionally written in terms of the Euler angles (ref, 17). However, since we are 
dealing with rigid rods, only two angles are needed to specify the orientation and these 
might just as well be chosen as the spherical angles 0 and <p (where e is the polar 
angle and (p the azimuthal angle). Hence, 


rot 


a).ij 


0? + sin^e. 




i = i,2 


(19) 


This form possesses the distinct disadvantage that it leads to the appearance of sin 0. 
in the denominator of the differential equation determinii^ (py resulting in imdesirable 
numerical consequences when 0j is near either zero or tt. This difficulty can be easily 
circumvented by describing the rotational motion of a rigid rod in terms of Cartesian 
coordinates 5^ (where j = 1,2) and by constraining these coordinates to lie on a unit 
sphere (ref. 18). The angles 0j and cp^ now describe the orientation of the unit 
vectors ^ . 

The foregoing considerations lead to 

(20) 

subject to the constraints 

^j^^j^ = i ^ Xj - 1) = 0 j = 1, 2 (21) 

The factor 1/2 in equation (21) was inserted simply for later convenience. If 
R = (X,Y,Z) is the separation vector for the two molecules and m their reduced mass, 
the relative translational kinetic energy is 


^trans 


R • R 


while the dipole-dipole interaction has the form 


( 22 ) 
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In this section juj and jH 2 are the molecular dipole moments and do not designate mo- 
ments of a distribution. Note that we have assumed that the dipole moment lies along the 
rigid rod. To incorporate the constraints it is convenient to define a new Lagrangian 
containing imdetermined multipliers Aj (where j = 1,2) by 

if S L + + X2?72 


and to consider the variational principle 
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if dt = 0 


(25) 


In the variation of this integral, R, x., and are treated as independent and the varia- 
tions in R and are required to vanish at the end points. This is just a special case 
of the problem of Lagrange in the calculus of variations (refs. 19 and 20). Carrying out 
the required variations gives 



Using the fact that 


5R = -i 6R 
dt 


and performing an integration by parts , 




The integrated terms disappear because the variations vanish at the end points. Thus the 
variations in R and lead to the equations 



while the variations in give the constraint relations (eq. (21)). 

The rather simple structxire of the dipole-dipole interaction and the constraint rela- 
tions makes it possible to avoid the explicit appearance of the multipliers in the equations 
of motion. Since 9 t7j/ 03^ = ^, it follows from the constraint relations (eq. (21)) that 
Xj • (0 tjj/ 9^) = 1. Equivalently, this follows from Euler’s theorem on homogeneous func- 
tions. Hence, from the second member of equation (26) 



12 



where again Euler’s theorem on homogeneous functions was used. Also since 
3L/3^ = and by twice differentiating the constraint equations (21) with respect to 
time to obtain • ^ + Xj • = 0, we find that 

dt \9x.y 

Substituting equations (28) and (29) into equation (27) yields 

^ = V - 2Tr„t(j) 


( 29 ) 


(30) 


The complete set of working equations is now obtained by using equation (30) in equation 
(26), and adjoining equation (21) and its first derivative to the system to obtain 


MR = 

- 




0R 

rot^^^l 


av 



ax. 

^ • 


= 1 

V 


= 0 


i = 1,2 




(31) 


J 


These equations are only applicable if (R • R) > ct ; when (R • R) = o- , the molecules 
tindergo a hard -sphere collision. The derivatives in equation (31) have the specific forms 

|(xj* Xg - 5r • XjX 2 * + r • Xj^Xg + XjXg- rl 

0R (R- R)2 ^ ^ 


„ ) 

ax, (i- 5)3/2 


av ^1^2 


^2 (R‘ 


- 3lf?-Xj) 


(32) 


where r = R/(R‘ R)^^^. 
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The equations of motion (eq. (31)) are solved numerically with a variable-step-size, 
fifth-order Runge-Kutta integration scheme. The method of altering the step size is de- 
scribed in appendix A. At each step of the integration the component of ^ with the 
largest absolute value is determined. Its value and the value of its time derivative at the 
end of the next step are calculated from the last two members of equation (31); concom- 
itantly, the differential equation for this component is ignored. The remaining variables 
and their time derivatives are obtained from the differential equations contained in equa- 
tion (31). The change in step size was determined by an empirical technique based on 
the error accumulated in the total energy of the pair of molecules. In all cases the total 

4 

energy is conserved to better than about one part in 10 . Random checks of time rever- 
sal showed that the coordinates returned to within several tenths of a percent of their ini- 
tial values. The velocities were within a few percent of the initial values except for very 
few trajectories. The integration was terminated when the molecules were sufficiently 
separated so that the interaction potential was effectively negligible; that is, 

< min [Tj.q^( 1), T^.^^ (2)j. The average calculating time on an IBM 7094 
was about 8 seconds per collision, although those trajectories characterized by low rela- 
tive translational energies required considerably more computer time. 

Each set of initial conditions for a potential collision was examined to determine 
whether the trajectory should be calculated. There were three circumstances under 
which a set of initial conditions was rejected. The first category of rejected initial con- 
ditions is that corresponding to collisions with projected running times in excess of 1 
minute. A second class of rejected initial conditions is that corresponding to collisions 
with large impact parameters b. Such collisions can reasonably be expected to produce 
negligible changes in the rotational energy. These collisions were rejected if b > 2. 5 a, 
where a was the same for all molecules; this corresponds to impact parameters 
b > 8.42 A. The third category of rejected initial conditions differs from the first two. 
The calculation of the trajectories in this category was attempted but had to be prema- 
turely terminated because the integration step size became too small. The minimum al- 

20 fi 

lowable step size used was 10“ seconds = 10“ scaled units. Collisions that fell into 
this category were generally those involving molecules with high rotational energies. 

Very little computing time is used to recognize these cases since the step-size changing 
routine rapidly reaches this minimum step size for these collisions. 

It is shown in the following section that it is convenient to give the initial rotational 
conditions of a molecule in terms of the angles 6 and 0 rather than the Cartesian coor- 
dinates ? (suppressing the index j). These quantities are related by the expression 
X = (sin e cos 4>, sin e sin (p, cos e). The time derivatives are then obtained by 
differentiation 

X = (0 cos e cos 0-0 sin e sin 0, e cos 6 sin 0+0 sin e cos 0, -6 sin e) 
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SELECTION OF INITIAL CONDITIONS 


Since the scattering problem is a two-body calculation, it can be reduced to a prob- 
lem involving only the rotational motion of the two molecules and their relative transla- 
tional motion. The pertinent equations are described in the previous section. A com- 
plete set of initial conditions for these equations requires the selection of a ’’sample’* 
molecule and a specification of its initial rotational motion, a specification of the initial 
rotational motion of a ’’heat bath” molecule, and a specification of their relative transla- 
tional motion. All these assignments of initial conditions must be made in a manner that 
is consistent with the fact that ( 1 ) ”heat bath” molecules are in an equilibrium distribu- 
tion at a temperatvire T, (2) the ’’sample” and ’’heat bath” molecules are in translational 
equilibrium, and (3) the ’’sample” molecule has a known rotational energy. 

First, the problem of the specification of the initial rotational state of a rigid -rod 
molecule is examined. If the effect of interactions is neglected, the rotational coordi- 
nates e and <p and their conjugate momenta (see eq. (19)) 


0T. 


rot 




de 


= Pg = le 


0T. 


rot „ c^v.2 n 

= = 10 Sin e 


00 


(33) 


have the normalized distribution 


‘^Pe «^P 0 = -^ exp 


8Ijr‘ 




21 \ ._2 


sin e 


de d0 dpg dp^ (34) 


where O:S0<it, Os 0 < 2t 7, -«> ;< < oo^ _oo ^ p^ :<oo. if we introduce the 

tr ans formation 


p^ = •^ p sin e cos 0 




p P sin 0 




(35) 


e = 0 


0=0 
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with Jacobian Ip sin e, the distribution takes the simpler factorized form 


1-^ h> 512^6^ (36) 

I "I \ ^ / ^ 

for 0:sp:^°o, 0 < < 2v, Q ^ 6 and 0 ^ 4> ^ 2n. The new variable p effectively 

specifies the rotational energy, while relates the two momenta p^ and p^. Thus, 
from the first two members of equation (35), we readily deduce 



.1 /Pfl ^ 
»// = tan ^ 1-2 

\ P0 ^ 


Each factor on the right is now a normalized one-variable distribution function for which 
the cumulative distribution (ref, 15, p. 23) can readily be obtained by integration. 
Therefore, we can pick at random from this distribution by using the method described 
by Schreider (ref. 21). In this method the cumulative distribution is equated to a random 
number R (numbers which are uniformly distributed over the interval [0, 1]), and the 
equation is then solved for the independent variable. This leads to 


-(f) 


In Rp 


x{/ = 2 ti R, 




e = cos"^ (1 - 2Rg) 


<p=2nR^ 


(37) 


where in the first of these equations we used the fact that if R is a random number, so 
too is 1 - R. The subscripts on the various R’s serve only to distinguish different ran- 
dom numbers. By combining equation (33) with equations (35) and (37), expressions are 
obtained for the quantities 0 and (p which are used directly rather than p^ and . 
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( 38 ) 



Equation (38) and the last two members of equation (37) can be used to assign the ro 
tational state of a randomly selected ’’heat bath” molecule. The orientation of a ’’sam- 
ple” molecule may still be calculated from the last two members of equation (37); how- 
ever, equation (38) cannot be used for their velocities since the rotational energy of 
the i*^ ’’sample” molecule is known. Nevertheless, since p^ is just twice the rota- 
tional energy, we may replace equation (38) with 




1/2 


e, = 


sin 27 tR, 





2e, 


nl/2 


IR,(1 - V 


> 

cos 2ttR^ 


(39) 


for the ’’sample” molecules. 

There still remains the task of assigning the initial coordinates and momenta for the 
relative translational motion. The ’’sample” molecule is taken as the origin of the Car- 
tesian, orthogonal coordinate system describing the relative motion (see fig, 3), while 
the ’’heat bath” molecule is situated on the negative Z-axis at a distance p from the 

o 

origin. The initial separation distance was chosen as p = 10 A. Let Pj, P 2 , and Pg be 
the translational momenta conjugate to X, Y, and Z, respectively. Obviously, if any 
scattering is to take place, pg must be greater than zero. Thus, it is appropriate to 
pick the p. from the normalized flux distribution 


^trans^Pl’P2’P3>‘^Pl ^P2 *^P3 



exp 


-^(Pl + P2 + Pg) 


2m 


Pg dp^ dpg dpg 


(40) 


for 0 Pg < 00 , -oo < :< oo^ and -°o < pg ^ If the transformation 
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for O^q^oo, O^p’^oo, and 0 ^ <p ^2 it. The variable q is just that part of the di- 
mensionless relative kinetic energy that arises from the velocity components in the X-Y 
plane. The variable ]p is that part of the dimensionless relative kinetic energy that 
comes from the Z -component of velocity 

2m 

^ ~ ^^trans " 

Once again, the cumulative distribution function method can be used for picking from 
^trans- This gives 

q = (- In 
p= (- lnRp)l/2 
<P=2nR^ 

where in the first two equations we again used the fact that if R is a random number, 
then so too is 1 - R. Combining equations (43) and (41) gives 
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Pl = 
P2 = 



( 44 ) 


This completes the specification of the initial conditions for the scattering problem apart 
from the selection of the particular ’’sample" molecule that is to participate in the scat- 

■Hi 

tering. If the r“ molecule is to participate in the scattering, the index i is chosen by 
the prescription 


"1 + [N%] % < 1 

i =< 

= 1 

where [a] represents the largest integer less than a. 


(45) 


DESCRIPTION OF COLLISON TRAJECTORIES 

Time-history plots of relative translational velocity and target rotational energy and 
collision trajectories are shown in figures 4 and 5 for two types of collisions. The plots 
in figure 4 correspond to a reflecting collision, that is, a collision for which the distance 
of closest approach is a. The plots in figure 5 correspond to a simple deflection colli- 
sion, that is, a collision for which the distance of closest approach is greater than a. 

All these plots were constructed by connecting N + 1 points, corresponding to N inte- 
gration steps, by N straight lines. The plotting was terminated after no more than 250 
integration steps. 

A typical variation of the absolute value of the relative velocity is shown in figure 
4(a) for a reflecting collision with an initial value of |^| nearly thermal at T = 300 K. 
The corresponding variation in rotational energy of the target molecule is shown versus 
intermolecular separation in figure 4(b). The periodic character of both plots is due to 
the change in orientation angle between the dipoles. This oscillatory behavior becomes 
more irreg^ular and of higher amplitude at small separations, where the dipole-dipole 
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interaction becomes strongest. Since the relative translational energy mv /2 = 0. 016 eV 
is much greater than the rotational target energy e « IcTq = 0. 002 eV for these condi- 
tions = 300 K; = 25 K), the change in |r| (~50 percent) is much smaller 

than the change in € (a factor of 8). The projections of the incident heat bath molecule in 
the three Cartesian planes are shown in figures 4(c) to (e) for the simple reflection colli- 
sion where the radial velocity changes sign at 3. 36 A. Similar plots are shown for a sim- 
ple deflection collision in figures 5(a) to (c). The variation in initial velocity is relatively 
small since the turning point for relative motion is at 4. 6 A. The target rotor is still 
considerably perturbed since it is initially relatively ’’cold”; the rotational energy varies 
by about a factor of 5. The projection of the incident molecule’s trajectory in the Y-Z 
plane is nearly a straight line. The projection of the motion onto the other two coordi- 
nate planes also gave straight lines. 


RESULTS 

Calculations have been performed for both HCl and DCl at two temperatures, 300 
and 500 K, These foiar calculations were supplemented by 300 K calculations for two 
pseudo -HCl molecules which had the same molecular properties as HCl except for the 
dipole moment. In one case, the dipole moment was one-half the dipole moment of HCl; 
while in the other case, the dipole moment was 9/10 the value for HCl. These two 
pseudo-HCl molecules are designated HCl(l/2) and HC1(9/10). The molecular properties 
for the four molecules are listed in table I. Each of the six calculations was carried out 
for 3000 collisions on a sample of 50 molecules drawn at random from an equilibrimn 
distribution of rotational energy at Tq = 25 K. The same initial distribution was used in 
all six cases and corresponded to a mean initial rotational energy of /3 qEq = 0. 9464. 

Equation (4) was fitted to the results of the calculations both as a two parameter 
equation and a one parameter equation by the methods described in the appendix B. The 
two parameter fit was always superior to the one parameter fit, just as would be ex- 
pected. The fitting was carried out after every 50 collisions, and the results are pre- 
sented in figures 6 to 11. It is apparent that the choice n^ = 2000 will satisfy the con- 
dition nj > for all cases except for HCl(l/2) in the one parameter fitting of fig- 

xire 11. In spite of this, n^ = 2000 is used even in this case. Table n presents the mean 
values of the fitting constants in the interval [2000, 3000], as well as the initial relative 
kinetic energy and the initial rotational energy of the heat bath molecules E^^^, 

averaged over the 3000 collisions. Five of the six values of /3E^ given in table H lie 
within 20 percent of unity, which agrees astonishingly well with the prediction of 84 per- 
cent given in figure 2. Similarly, the values of are all within 5 percent of the 
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infinite population average of 2. The fact that is consistently greater than the 

infinite population average of 2 is most likely a reflection of bias introduced by the proc- 
ess of rejecting initial conditions. The most obvious culprit of the three rejection cat- 
egories is the rejection based on excessively long projected running times. This con- 
sistently rejects only collisions with low relative kinetic energies. This category repre- 
sented 20 to 30 rejects for every 3000 collisions used. It is unlikely that any bias was in- 
troduced by the large-impact -parameter rejection category since this criterion only ex- 
presses a relation between the component of K along the Z -axis and the other two com- 
ponents. Since the components are assigned randomly, this should reject some low rela- 
tive kinetic energy collisions as well as some high energy collisions. However, initial 
conditions rejected because the step size was too small undoubtedly played a part. This 
was particularly obvious when we attempted to do calculations for HF. We found that, 
because of the strong dipole-dipole interaction, approximately 20 to 25 percent of the 
collisions that we attempted to calculate were being rejected becaiise the step size be- 
came too small. This produced such a strong bias in that we considered the 

calculations imusable. Typically, there were 800 to 1000 large-impact-parameter re- 
jects, 20 to 30 long- integration-time rejects, and 20 to 30 step-size-too -small rejects 
for the calculations reported here. 

Figures 12 to 14 display two exponential curves calculated from equation (4), as well 
as their asymptotic values and the data points for Epj after every 10 collisions. One 
curve was obtained using /3E^ = 1 and the mean one parameter for the interval 

[2000, 3000]. The other curve corresponds to the mean two parameter for the in- 

terval 2000 to 3000 and the /3E^ which minimizes the function # of appendix B for this 
Zrof Table in tabulates the function for each pair of curves and also the value of # 
obtained by a direct one or two parameter least squares fittii^ of the 3000 data points. 

The data in table O clearly demonstrate the superiority of a two parameter fit over a 
one parameter fit since the value of 4> is always smaller for a two parameter fit. They 
also show that the mean Z^^^. curves are only moderately inferior to the free Z^^^ 
curves. A cursory examination of figures 12 to 14 discloses that an exponential does 
indeed seem to describe the evolution of the sample rotational energy but that the data 
scatter considerably about the fitted curve. 

Table rv contains a summary of the experimental and theoretical for HCl and 

DCl. The experimental numbers come from both thermal conductivity measurements and 
acoustic measurements. Even though the validity of the thermal conductivity numbers 
was questioned in reference 2, we nevertheless include them for comparison purposes. 
The acoustic numbers are the recent results obtained by Evans (ref. 22). He analyzed 
his measurements in terms of both a single and a double relaxation process and, further, 
he used a frequency dependent specific heat when calculating the classical absorption. 

The theoretical Z^^^ were obtained from tabulations in reference 2, as well as from the 
calculations of this report. In the molecular dynamics results we have included not only 
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the mean for the interval [2000, 3000], but also the extreme values in this interval. 

The latter are a more meaningful measure of the uncertainty in this case than the stand- 
ard deviation. 

Several qualitative observations can be made on the basis of the data in table IV. 
First, all the molecular dynamics results fall between the experimental data and the re- 
sults of the perturbation calculation and, hence, are in closer agreement with experi- 
mental results. Second, all the numbers in table IV show that Z^^^ increases with tem- 
perature except for those extracted from thermal conductivity measurements and the one 
parameter results for DCl. This is also apparent from table V. Third, all the theoret- 
ical numbers, except the one parameter, 300 K molecular dynamics data, show that in- 
creasing the moment of inertia decreases Z^^^. On the basis of these facts it seems 
that the two parameter data should be preferred over the one parameter results. 

Fourth, a comparison of HCl and HCl(l/2) shows that the dipole moment dependence of 
the molecular dynamics results is considerably weaker than the dependence obtained 
in references 1 and 2; this is particularly true of the two parameter numbers. 

We have shown that the results from the two parameter fit of the molecular dynamics 
calculations are in reasonable consonance with other theoretical predictions and, in addi- 
tion, seem to be a closer approximation to the acoustic data than were the perturbation 
calculations of references 1 and 2. But this alone is Insufficient evidence to inspire an 
overwhelming confidence in the results. Any such confidence is quickly tempered by the 
incongruous result for HC1(9/10) in table IV. Based on the values for HCl and HCl(l/2) 
we would expect Z^^^ for HC1(9/10) to be quite close to that for HCl but still slightly 
larger. It is in fact considerably smaller. This opens two possibilities for considera- 
tion. The first is that 3000 collisions are insufficient to obtain an accurate Z^^^; the 
second possibility is that a 50-molecule sample is too small to give a precise value for 
Zrof The first option can be easily investigated by determining how close the fitted 
curves for Epj, shown in figures 12 and 14, are to their asymptotic values after 3000 col- 
lisions. This is shown in table VI. Obviously, in all cases 3000 collisions brings the 
sample within 91 percent of the asymptotic value for Epj except for the one parameter 
data for HCl(l/2). As a matter of fact, HCl and HC1(9/10) at 300 K are each within 5 per- 
cent of their asymptotic value. This certainly must eliminate the first possibility and 
leave us only with the second one. So it appears that even under the best of circum- 
stances we must be prepared to accept an uncertainty in Z^^^ of the order of six units, 
representing the difference between HCl and HC1(9/10) when working with such a small 
sample. 

Finally, there is one additional aspect of the problem of determining Z^^^ by molec- 
ular dynamics. To carry out the scattering calculations we were required to select an 
initial separation p. Our choice of an initial separation of p = 10 A was dictated by 
pragmatic considerations. However, this choice committed us to the proposition that any 
potential collision with an impact parameter b > p is not a physically significant 
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collision insofar as gas phase properties are concerned. Thus, in effect, we made a de- 
cision on what should be called a collision. However, in our computations we did more 
than this, since we rejected all collisions with impact parameters that exceeded 
bm = 2. 5 a = 8.42 A. Thus, collisions whose impact parameters were in the range 
8.42 A < b < 10 A were not counted. Although we cannot estimate the effect of our 
choice for p, we can quite easily place an upper limit on the effect of our choice for b^^. 
To do this we must know what fraction s of all physically significant collisions (b < p) 
have impact parameters less than, or equal to, b^^. But this is nothing more than the 
cumulative distribution function for b evaluated at b_ . 

The impact parameter is related to the variables q and p of the translational dis- 
tribution function (42) by 


/u\2 2 

/ b r _ q 

V n) 2 —2 

q +P 

2 —2 

Thus, in the space spanned by q , p , the curves of constant impact parameter are the 
straight lines 





Thus , the cumulative distribution function for b is given by the integral 


I d(p2) I 

•'0 Jo 


2 2 

Therefore, it follows that s = b^/p . If we now reasonably assume that collisions with 
b >bj^ are uniformly distributed during the course of the calculation, the choice of 
bvv, < p implies a uniform scale change in the abscissa n, n = n's. Hence, since 
n/ NZrot = we conclude that = p‘^ Incorporating our 

choice of p and b^^^ results in =1.41 Z^^^. This is clearly an upper limit to the 

effect since this presumes that the ordinate Ej^ is unaffected. In spite of the fact that 
the neglected collisions would not be responsible for a large fraction of the total energy 
transferred to the sample molecules, they could transfer some energy. This is partic- 
ularly true during the initial phases of the calculation, when the net effect should be to in- 
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crease the energy absorption by the cold sample molecules. This could somewhat com- 
pensate for the stretching of the abscissa. Interestingly enough, the effect we have been 
discussing serves to bring the rotational collision numbers somewhat closer to the values 
obtained by the perturbation calculation. This is shown in table VII which gives 
Z^ot» numbers obtained by the perturbation calculation for HCl and DCl. Al- 

though, the values of Z^^^ are closer to the perturbation results, they are further away 
from the experiment acoustic values than ^rof 

CONCLUSIONS 

The discussions of the previous section can be summarized as follows: 

1. The mean rotational energy of the sample molecules exhibits an exponential ap- 
proach to a steady state value; however, there is considerable scatter about the expo- 
nential curve. 

2. The molecular dynamics results for HCl and DCl are generally in consonance with 
previous perturbation calculations. 

3. The values of for HCl and DCl obtained by molecular dynamics are 

intermediate to the experimental Z^^^ values and those obtained by a perturbation 
calculation. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 10, 1971, 

129-01. 
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APPENDIX A 


STEP-SIZE CONTROL 


The integration step size was determined by an empirical technique based on (1) the 
conservation of energy E during a collision, (2) a specified maximum allowable error 
in the energy, AE > 0, and (3) an estimated total trajectory time T. If E^®^ > 0 is the 
initial energy, h^^^ is the step size to be taken in the next integration step, 

= t^^^ + h^^\ and E^^^ and are the values of the energy before and after the 

step, respectively, then the step produces an actual error increment in the energy that 
may be expressed as 



^( 0 ) 


On the other hand, a maximum tolerated error increment for the step is defined by 


(Al) 


g(i+l) ^ [ aE - |E^^) - E^°) h£_^ 

e(o)t 


(A2) 


At the conclusion of each step we calculate 


^(i+1) ^ ^ - E^^^ I 

AE - |E^^^ - E^®) I 

0 < < 2, the step is accepted; otherwise, it is rejected. In either case, a new 

step size is chosen by the prescription 


T 

h(i) 


(A3) 





(A4) 


where m is the order of the Runge-Kutta integration routine. 

This prescription for altering the step size can lead to either a larger or a smaller 
step size. Should the step size become smaller than 10"^® second, the calculation is 
terminated. Fortunately, such a situation rarely occurs. 
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APPENDIX B 


LEAST SQUARES FITTING 

There are three available parameters to fit equation (4) to the sample rotational 
energy data in a least squares sense, namely, Eq, E^, and We chose Eq so 

that the curve would pass through the known rotational energy prior to any collisions. 
This leaves only E^ and The parameter E^ presumably corresponds to the 

equilibrium value which is equal to kT for an infinite population of molecules with two 
rotational degrees of freedom. But for a sample of only 50 molecules, based on our sta- 
tistical considerations (fig. 2), there is no assurance that E^ will turn out to be kT. 
Thus we determined both E^ and by the least squares fitting procedure. 

At this point it becomes convenient to introduce some alternate notation. Thus, we 
define 


g = /3(Ej^ - Eq) 


and introduce new parameters by 


A = /3(E^ - Eq) 




X = Z 


-1 

rot 


Then equation (4) becomes the function 


^ (n) = A 


1 - exp 


(-3 


If we denote the tabulated data by ^(n), we can define 


4-(x,A)=l [,r(n) - g{n)f 

2 


n=] 


(Bl) 


(B2) 


(B3) 


(B4) 


where M is the total number of points used in the fitting. The expression for the sum 
of a geometric series can be used to reexpress in the form 


26 



#(x, A) = [M + 1 - 2g(x) + g(2x)] - A 
2 


M 

J(n) - J(n)e~ 

n=l n=l 


nx/N 


M 


+ - (B5) 


n=l 


where 


1 - exp 


g(x) = 


r (M + i)x 

L N ^ 


--[•I] 

The function g(x) is monotonically decreasing and has the limiting values 

g(o) = M + r 

g(oo) =1 J 

Thus, it follows that # has the limiting values 


M 

^(0,A) =1 ^ J^{n) 


M 




#(oo, A) = 1 MA^ - A ^(n) +~y^ J^(n) 

2 ^ 2 

n^ n=] 


(B6) 


(B7) 


(B8) 


The parameters A and x are determined so as to minimize that is, they are 
solutions of the equations 
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— = A[M + 1 - 2g(x) + g( 
0A 


M , 

(2x)]-^ 5(n)(l - = 0 


M 

I 

M 


Na“ 1 ^ = AN[g’(2x) - g’ (x)] - \ ^ n J(n) 
dx 


7m «-nx/N ^ 


= 0 


n=l 


(B9) 


where the primes are used to symbolize differentiation with respect to x. The first of 
these immediately gives A as a function of x. 


A(x) = 


M 

E J(n)(l - 

n=l 


M + 1 - 2g(x) + g(2x) 


(BIO) 


We now define a function of x alone as 


G(x) = NA"^ — 
dx 


A=A(x) 


(BID 


The solution of the equation 


G(x) = 0 


(B12) 


can now be accomplished by the Newton -Raphson iteration (ref. 23) 






(B13) 


where 0 < X < 1 is a convergence control parameter. The iteration is carried out by 
using the following expressions: 
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M 


g’(x) = A’Njg'(2x) - g’(x)j + ANj2g"(2x) - g"(x)j + n^(g’(n) 


„-nx/N 


n=l 


N 


A’(x) = =* 


_ -nx/N 

n,?(n) + 2[g (x) - g (2x)]A 

N 


M + 1 - 2g(x) + g(2x) 


l.e-VN 


Ng"(x) 


M+ 1 _ fM H. 1 - 

1 - (, . e-VN)2 


The actual procedure used was to make an Initial estimate for x, calculate A(x) from 
equation (BIO) and obtain a new estimate for x from equation (B 13). The iteration was 
terminated whenever | Ax^^^ | < 10”®. 

When was not to be treated as a fitting parameter, we set A = 1 - (EQ/kT) and 
A’ = 0 and proceeded as before. When was not to be treated as a fitting parameter, 

A could be calculated directly from equation (BIO). 
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TABLE L - MOLECULAR PROPERTIES 


[Molecular diameter^, 3.36x10"^ cm (3.36x10"^^ m). ] 


Molecule 

Molecular 

weight 

Moment of 
inertia, 

I, 

g-cm 

Dipole 

moment, 

At, 

esu-cm 

HCl 

36.461 

2.6431x10“^^ 

1.081x10“^® 

DCl 

37.467 

5. 1404 

1.085 

HCl(l/2) 

36.461 

2.6431 

.541 

HC1(9/10) 

36.461 

2.6431 

.973 


^All molecular diameters set equal to that of HCl, which 
was obtained from an analysis of viscosity data made 
by using the Stockmayer potential (ref. 24). 


TABLE II. - MEAN FITTING PARAMETERS IN INTERVAL 2000 TO 3000 AND MEAN RELATIVE KINETIC 
ENERGY AND HEAT BATH MOLECULE ROTATIONAL ENERGY FOR 3000 COLLISIONS 


Molecule 

Temper- 

Two parameter fit | 


ature, 

Dimensionless 

Rotational 


K 

rotational 

collision 



energy, 

number. 




^rot 

HCl 

300 

1. 17 

19.0 

HCl 

500 

1,07 

24.4 

DCl 

300 

.97 

14.2 

DCl 

500 

1. 19 

21.7 

HC1(1,'2) 

300 

.64 

22.3 

HCl (9 10) 

300 

1.02 

12.9 


Rotational 

Mean relative 

Heat bath 

collision 

kinetic energy, 

molecule 

number 

^rot 

one parameter 
fit 

^^trans 

rotational 
energy for 
3000 collisions, 

^^rot 

12.9 

2. 10 

0.98 

21. 1 

2.07 

.98 

15.6 

2.06 

1.00 

14.4 

2.03 

.99 

51.7 

2.04 

.97 

12.2 

2.08 

1.00 


TABLE III. - THE FUNCTION 4> AFTER 3000 COLLISIONS 


Molecule 

Temper- 

ature, 

K 

# for two parameter fit, 
rotational collision 
number - 

$ for one parameter fit, 
rotational collision 
number 

Free 

Fixed at mean 

Free 

Fixed at mean 

HCl 

300 

10.462 

10.478 

18.093 

18. 138 

HCl 

500 

2.883 

2.902 

3.667 

3.703 

DCl 

300 

7.083 

7.357 

8.727 

8.774 

DCl 

500 

4.458 

4.910 

17.374 

17.537 

HCl(l/2) 

300 

1.236 

1.446 

2.576 

2.693 

HCl (9/ 10) 

300 

10.983 

12.256 

14.747 

14 . 803 
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TABLE IV. - EXPERIMENTAL AND THEORETICAL ROTATIONAL COLLISION 


NUMBERS FOR HCl AND DCl 


(a) Experimental 


Molecule 

Temperature, 

Rotational collision number Z^^^ 

obtained from- 


K 





Thermal 

Acoustic measurements^ 



conductivity 

measurements^ 

Single - 
relaxation 

Double- 

relaxation 




process 

process 

HCl 

300 

3.8 

6.4 

3.0 

HCl 

500 

1.0 

6 

‘^4.3 

DCl 

300 

.9 j 

— 

— 

DCl 

500 

0 

— 

— 

HCl(l/2) 

300 



--- 

HCl (9/ 10) 

300 

--- 


— 


(b) Theoretical 


Molecule 

Temperature, 

K 

Rotational collision number 

calculated by- 



Perturbation 


Molecular dynamics® 





calculation^ 

Two parameter fit 

One parameter fit 




Mean 

Max. 

Min. 

Mean 

Max. 

Min. 

HCl 

300 

65.9 

19.0 

21.2 

17.5 

12.9 

13.2 

12.6 

HCl 

500 

79.8 

24.4 

27. 1 

23.2 

21. 1 

21.4 

20.8 

DCl 

300 

51.9 

14.2 

16. 1 

12.7 

15.6 

16.0 

15.3 

DCl 

500 j 

64.2 

21.7 

25.0 

18.7 

14.4 

15.0 

13.8 

HCl(l/2) 

300 

— 

22.3 

27.4 

19. 1 

51.7 

53.4 

48. 7 

HC1(9/10) 

300 i 

— 

12.9 

16.4 

10.0 

12.2 

12.5 

11.8 


^Taken from table ITT of ref. 2. 
^Data of L. B. Evans (ref. 22). 


^Read from a curve faired through data of L. B. Evans (ref. 22). 
^Taken from table IV of ref. 2. 

^This work. 
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TABLE V. - TEMPERATURE DEPENDENCE OF 


ROTATIONAL COLLISION NUMBER 
FOR HCI AND DCI 


Source 

Zj.^^(300 K)/Zj.^^(500 K) 


HCI 

DCI 

Experimental: 



Thermal conductivity 

3.8 

DO 

Acoustics 



Single -relaxation process 

,67 

— 

Double -relaxation process 

.70 

— 

Theoretical: 



Perturbation calculation 

.83 

0.81 

Molecular dynamics 



Two parameter fit 

,78 

.66 

One parameter fit 

.61 

1.08 


TABLE VI. - APPROACH OF SAMPLE ROTATIONAL ENERGY Ej^ 


TO ITS ASYMPTOTIC VALUE AFTER 3000 COLLISIONS 


Molecule 

— 

Temperature, 

K 

CEj^(3000) - Eg] (E„ - 

Two parameter fit 

One parameter fit 

HCI 

300 

0.9572 

0.9903 

HCI 

500 

.9142 

.9415 

DCI 

300 

.9855 

.9787 

DCI 

500 

.9372 

.9845 

HC1(1 2) 

300 

.9326 

.6865 

HC1(9 10) 

300 

.9904 

.9927 


^Calculated by eq. (4) using the data given in figs. 12 to 14. 
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TABLE VII. - COMPARISON OF THEORETICAL 


ROTATIONAL COLLISION NUMBERS 


Molecule 

Temperature, 

K 

Rotational collision number, determined by- 


Perturbation 

Molecular dynamics 



calculation 

Two parameter 

One parameter 




fit 

fit 




^rot 

^rot 

^rot 

Kot 

HCl 

300 

65.9 

19.0 

26. 8 

12.9 

18.2 

HCl 

500 

79.8 

24.4 

34.4 

21. 1 

29.8 

DCl 

300 

51.9 

14.2 

20.0 

15.6 

22.0 

DCl 

500 

64.2 

21.7 

30.6 

14.4 

20.3 


Exponential approximation of ensemble evolution 



a> 



Figure 1. - Temporal evolution of rotational energy of ensemble. 
Shading represents density of data points. 





.-\Dsolute ratio of velocity to initial velocity 




Absolute ratio of velocity to initial velocity 


1.015 

1.010 

1.005 

1.000 

.995 

.990 

.985 

.980 

.975 




(c) Collision trajectory in Y-Z plane. 

Figure 5. - Time-history plots of relative translational velocity, target rotational energy, and collision 
trajectory for a simple deflection collision. 








HCI(l/2)-300 K 



loloo ^ 1500 ^ 20*00 

Number of collisions, n 


2500 


3000 


Figure 8. - Least squares fitting constants for HCHl/2) and HCK9/10) as a function of the number of collisions; 
two parameter analysis; temperature, 300 K. 
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Lire 9. - Least squares fitting constant for HCI as a function of the number of collisions: one parameter analysis. 
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Figure 10. - Least squares fitting constant for DCI as a function of the number of collisions- one parameter 
analysis. ' ^ 


3( 



HCI(9/10)-300K 



Number of collisions, n 


Figure 11. - Least squares fitting constant for HCKl/2) and HCK9/10) as a function of the number of collisions: 
one parameter analysis. 
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(a) HCK9/10). 
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(b) HCKl/2). 

Figure 14. - Rotational relaxation of HCK9/10) and HCKl/2). Temperature, 300 K. 
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